Iranian Journal of Numerical Analysis and Optimization 
Vol. 14, No. 1, 2024, pp 136-171 


https: //doi.org/10.22067 /ijnao.2023.82954.1281 YU 


https://ijnao.um.ac.ir/ 


Research Article 
BY NC 


A study of global dynamics and 
sensitivity analysis of a discrete-time 
model of the COVID-19 epidemic 


M. E-B. Keddar and O. Belhamiti*® 


Abstract 


This study presents a novel approach to understanding the global dynam- 
ics of COVID-19 transmission with vaccination based on a discrete-time 
model. We establish biologically meaningful constraints on the model 
parameters and prove the existence of a disease-free equilibrium and an 
endemic equilibrium under these constraints, along with their theoretical 
stabilities. Furthermore, we identify the most sensitive operational param- 
eters that have a substantial impact on the transmission of the epidemic. 
Through numerical simulations, we demonstrate the local stability of the 
two equilibria, depending on the parameter values. Our findings reveal 


that the discrete-time model is not only dynamically robust but also more 
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realistic than its continuous counterpart under biologically meaningful con- 
straints. These results provide a foundation for future research in this area 
and contribute to our understanding of the global dynamics of COVID-19 


transmission. 


AMS subject classifications (2020): Primary 39A12; Secondary 92C60, 92D30. 


Keywords: Discrete-time epidemic model; Euler Scheme, Local stability; 


Vaccination reproduction number. 


1 Introduction and Motivation 


Coronaviruses are a family of viruses that cause respiratory and intestinal in- 
fections in humans and animals [6]. They were not dangerous to humans until 
the appearance of severe acute respiratory syndrome coronavirus (SARS- 
CoV) in 2002 in China [16]. This virus had infected 8,096 people in 25 
countries, causing 774 deaths [19]. A decade later, the outbreak of another 
highly pathogenic coronavirus, Middle East respiratory syndrome coronavirus 
(MERS-CoV), occurred in Jeddah, Saudi Arabia [32]. The virus spread, in- 
fecting more than 2,400 people in 27 countries and causing over 850 deaths 
[13]. In December 2019, a third of these dangerous coronaviruses appeared in 
the city of Wuhan, China. Forty-four cases were reported in this country on 
January 3, 2020 [24]. The virus migrated to neighboring regions and coun- 
tries. In less than a month, 106 confirmed cases were detected in 19 different 
countries, and the number of infected individuals reached 9,826, with 213 
deaths. The World Health Organization (WHO) declared the outbreak to be 
a public health emergency of international concern [24]. The Severe Acute 
Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is the virus that caused 
this disease, which has spread throughout the world. It was named COVID- 
19 by the WHO [23], which declared that this malady can be characterized as 
a pandemic [25]. On April 4, 2020, the number of infected people exceeded 
1 million, with 56,986 deaths worldwide [26]. Many countries have gone into 
total confinement and closed their borders to reduce the spread of the virus. 


The number of infected people increases more and more rapidly over time. 
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In December of the following year, more than 75 million confirmed cases and 
1.6 million deaths have been reported globally [27]. As of November 2022, 
COVID-19 still exists, and its spread worldwide has resulted in 629 million 
confirmed cases and 6.5 million deaths [28]. The number is increasing even 


now. 


This pandemic has left the world feeling helpless due to its rapid spread 
and the lack of specific medicines or treatments available to combat it. As a 
result, many countries have implemented non-pharmaceutical measures, such 
as confinement, mask-wearing, social distancing, isolation, and treatment of 
infected persons, in an attempt to contain the virus. However, the effective- 
ness of these measures has proven to be limited. This has prompted epidemi- 
ologists to develop vaccines as a solution. Many of the developed vaccines are 
safe and effective in providing protection against severe forms of COVID-19, 
offering a glimmer of hope for an end to this pandemic that has impacted 
humanity for over two years. The COVID-19 pandemic has affected all sec- 
tors of health and security, including the involvement of mathematicians to 
control and stop its spread. Numerous mathematical models have been pro- 
posed, allowing for the control of the pandemic, prediction of its evolution, 
and estimation of key transmission parameters, such as the basic reproduction 
number Ro, as well as analyzing the impact of vaccines on disease dynamics. 
Most of these mathematical models are nonlinear continuous-time systems 
of the SEIR type, which study the dynamic transmission of COVID-19, such 
as those presented in [21, 12, 15]. Some models also include vaccination, as 


seen in [1, 31, 30]. However, only a few models are discretized-time models. 


It should be noted that discrete-time models present a great advantage 
compared to continuous-time models, especially for COVID-19, as the indi- 
vidual goes through several phases (Susceptible-Exposed-Infected ...) and 
the phases of the disease represent the step of discretization, which can vary 
between five days to a week, depending on the individual. Despite these con- 
siderations, only a few articles have been published in discrete time that take 
into account the disease phases. In [8], the authors proposed a SEIR, model 
with several characteristics. Firstly, the exposed subpopulation is infectious. 
They also incorporated a feedback vaccination control law on the suscepti- 


ble. Furthermore, the model is subject to delayed partial re-susceptibility, 
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meaning that a partial immunity loss in the recovered individuals occurs af- 
ter some time. By the same authors in [9], a new discrete SEIR model has 
been developed. This model considers the injection of two doses with dif- 
ferent delays and a possible average different effectiveness. It also considers 
that the exposed have a transmission rate exposed-susceptible, which may 
be eventually distinct from the infectious-susceptible transmission rate. In 
[2], the authors proposed a new SIR, model with vaccination. Their model 
has two equilibrium points, namely, a disease-free equilibrium (DFE) and an 
endemic equilibrium. By considering both the forward difference system and 
the backward difference system, some stability analyses of the DFE point are 


performed. 


This paper describes a new discrete SEIR model for COVID-19 trans- 
mission dynamics that include vaccination. The model assumes that only 
symptomatic and asymptomatic infected individuals can transmit the disease 
and that hospitalized patients are isolated. It also considers the possibility 
of reinfection with the virus, as the vaccine may not be fully protective. [4]. 
Some exposed and asymptomatic infected can be vaccinated; this is due to 
the inefficiency of the tests [29]. The vaccination reproduction number R, 
and the equilibrium points are calculated, and their sensibility and stability 
analysis are studied. Numerical simulations are performed considering vari- 
ous scenarios. We present detailed results in the form of graphs and tables, 


including new numerical results obtained for model validation. 


The paper is organized as follows. Section 2 introduces the mathematical 
model and provides a preliminary overview. In section 3, the mathematical 
analysis is carried out to investigate the stability of the proposed model. 
This section focuses on the study of equilibrium points and their stability 
properties. Section 4 presents a sensitivity analysis of the model to investigate 
how changes in the model parameters affect the system behavior. Finally, 
section 5 presents the results and discussions of the model simulations, which 


provide insights into the system’s behavior and its implications. 
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2 Mathematical model 


In this section, we use the explicit Euler scheme to obtain a discrete-time 
model for the transmission of COVID-19. Based on the epidemiological status 
and the timeline of COVID-19 infection in [18], the population is divided as 


follows: 
e S is the susceptible class, 
e Eis the exposed class, for example, the people in the incubation period, 
e V is the vaccinated class, 
e I* is the asymptomatic infected class, 


e T° is a class of individuals diagnosed with moderate COVID-19 symp- 
toms, which is further divided into two subclasses [*! refers to indi- 
viduals experiencing symptoms in their first phase (i.e., within a week 
of symptom onset), as described in [11, 17] while I*? refers to individ- 
uals experiencing symptoms in their second phase (i.e., roughly one 
week after the end of the first phase of symptoms), also documented in 
references [{11, 17]. That is, 


Io =1 41%. 


e I’ refers to the class of infected individuals who require hospitalization, 
including those who need intensive care. This class is further subdivided 
into three subclasses based on the interval since symptom onset: I! 
represents individuals who are hospitalized within a week of symptom 
onset, as detailed in references [11, 17], while I’? denotes those who 
are hospitalized around one week after the end of the first phase of 
symptoms. Lastly, I’ pertains to individuals who are hospitalized 
approximately one week after the end of the second phase of symptoms, 


also documented in references [11, 17]. That is, 


P= 1 +7 +1, 
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e R represents the class of individuals who have recovered from COVID- 
19. 


e D represents the class of individuals who have died from COVID-19. 


Week 1 Week 2 Week 3 Week 4 


r In do > Jz ds | Js 82 > | 
| | | \ J — 
uN uw | re 8 t 


Disease-free Incubation Post infection 
period period Infection period period 


Figure 1: Timeline of infection and states transition scheme description. 


By applying the explicit Euler scheme, we can obtain a discretized SEIR 
epidemic model that incorporates vaccination (see Figure 1). Considering 
the week-long interval required for individuals to transition between classes, 
the resulting model is as follows: 


S (t+ 1) =S(t) + uN — mS (0) + mV (0) + OR (1) — pS (t) 


as) OFM O +I?) 


(2 (4) FI (t) +1 (2) 


E(t+1) =E(¢) + a8 (t) a 
(61 + B2 + 63 +n3 +) E(t), 

V(t+1) = V(t) + mS (¢) + nsE (t) + nal* (t) — (n2 + yo + 2) V (2), 
I? (¢+1)=1°() + fE()-(mtn+HI@), 

Is. (¢+ 1) =I (t) + BE) —-OAi+¢e+y)I* (ft), 

1s (t + 1) = I? (t) + ArT (€) — (y2 + w) I? (4), 

11 (+1) =I" () + BsE(t) — (A2 +n) I" (2), 

Th2 (¢4+ 1) = 1" (¢) + Aod™ (¢) + gl! @) — As +73 +61 +H) 1 @), 
T's (¢ + 1) =I" (¢) + Ag” (t) — (ya + 62 + w) I" (2), 

R(t+1) =R(t) + V(t) + I? () + ya? (4) + 31"? (t) 

+4I"3 (t)— (0+ pH) RO), 

D(t+1) = D(t) 4 6,1" (t) + dgI"8 (t), 


(1) 
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We assume that the total population N remains constant, 
N(t+1)=N(t). 


Taking into account the previous assumption, we can normalize the variables 
in the model (1), 


sy =8 9, py =FO veQe, rw=*9 
mg=2 8 m@=EO, me-F 8, mye, 
my=EO, py=2O, p@=2. 
Because 
S(t)+E(#)+V(t)+ 1° (t) +I (#) 4+ I? (ft) 
+ 7 (¢) + 7" (¢) 4 78 (() + R(t) + D(t) = 1, 
then 


D(t)=1-(S(@+EL(t)4+V (41° (t)4+I° (0) 
+1° (t)+1™ (t) +1 (t) +1" (t) + R(t) - 


Then, we have the following normalized reduced system: 


S(t+1)=S(t)+u—mS (t) + nV (t) + OR (t) — WS (t) 
—aS (t) (1% (t) + I (t) + 1°? (t)), 
E(t+1) = E(t) + aS (t) (14 (t) + 1° (t) + 1° (t)) 
(81 + 62 + 63 +n3 +) E(t), 
V(t+1) =V (t) + mS (t) + ns (t) + mal® (t) — (n2 +70 +) V(t), 
I(t+1) =1¢(t{)4+ HE (t)— (mtu +p) 17 (t), 
PMe¢+l =P (t)+ PE t)-Aritethyl'(), (2) 
12 (t+ 1) = T°? (t) + ArT* (t) — (ga + we) L°? (8), 
(¢+1)=1" (0) + PEO) -QO2ty)I™(), 
rhe (f+1) = 172 (€) + AI" (£) + pl*! (t) — As +73 + 61 + w) I”? @t), 
TPs (t +1) = 18 (t) + Ag I”? (t) — (ya + 62 + ws) 18 (£), 
R(t +1) = R(t) + 0V (t) +114 (t) + yoT®2 (t) + y3I"? (t) 
+7418 (t) — (0+ 4) RQ), 
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where 
pL represents the natural mortality, 
a represents the transmission rate of COVID-19, 


(n:);-734 represent the vaccination rates for each class, 


(5;);-7-3 represent the rates at which exposed individuals become infected, 
(7i);-oa Tepresent the recovery rates from each class, 

(Ai);-73 Tepresent the rates of transition from one class to another, 
(5i)j=1,9 signify COVID-19 mortality rates per hospitalized group, 

8 and 72 denote COVID-19 reinfection rates, and 

y represents the hospitalization rate in the symptomatic phase J*'. 


It is essential to assume that all model parameters are positive. Further- 
more, some conditions on these parameters are necessary to ensure that the 
model has biological significance. In the following subsection, we present and 


illustrate these conditions. 


2.1 Positive invariant set 


The domain 


Cat(5 BV Ir 7 nye. 4 
O<S+E4+V4I*tI 477241" 41° +I8 +R <1} 


is a positive invariant set for system (2). 

To render the model biologically meaningful, it is imperative to exhibit 
the belonging of solutions from the system (2) into Q. This, in turn, hinges 
on ensuring that the combined outflow rates for each state do not surpass 1, 
as follows: 


meats, y2tu <1, 

Bi + Bo2+ P3t+na+u<l, A2+tu< i, 

nm +y+H <1, and 4 A3+73+61+" <1, 
mtyte<l, yatdo+u <1, 
M+e+u<l, O+py<cil. 


Let (Sn, En, Vn, 12, 181, 182, 1%, 122, 123, R,) € O. Then, we have 


MIO UIE OM POM DOM. FOR VOM 2 
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Snti = Sn +e—71Sn +72Vn + ORn — @Sn (1% + Int + 19?) — wSn 
>ut(—-(m+atp)) Sn > 0, 
Enti = En + aSn(I% + Int + In?) — (81 + B2 + 63 +73 + 1) En 


> (1 — (61 + 62 + 83 +73 + )) En > 0, 

Vn tmSn + 73En + nal? — (n2 +70 + 1) Vn 
> (1— (m2 +70 + B)) Vn 2 9, 

Tray =IR + Bi En — (nat + Ww) Ih 
>(1-(m+7+4+4)) 12 20, 

Ty =In t+ BeEn -Qiteth) ti 

>(1—- (Ari + ¢4+y)) In! > 0, 

Ty =I + AIn — 2 +) I? 
>(1-(y2 + 4)) In? S 0, 

Tht = IR! + B3En — (Ao tu) In 
>(1-(A2+p)) In? > 0, 

12, = In? + roTn? + pla} — (As +93 +61 +p) In? 
> (1-3 +43 + 61 + 4) In? > 0, 

IMS | = In + AgIn? — (ya + 62 + pw) In? 

> (1— (ya + 62 + p)) IR > 0, 
Rati = Rn t+70Va + W118 + y2lr? 
> (1- (0+ 1) Rn > 0. 


= 
| 


+ yg In? + yalh3 — (0+ 1) Rn 


Thus, 


Sn4it Engi t Vag +1044 + 1, +12 $1 +8 +I + Rng 2 0. 
On the other hand, the sum of the system equations (2) gives 


h h h: 
Snt1 + En4i + Vandi t+ Tigi + A +02 t+ Ih + a + a + Rost 


=pt+(l—H)(Snt+En+Vat+IP 


n 


7st 4 782 4 7h 4 he 4 P34 RR, 5, 1?2 + by 1h3 
nm n n n nr n n 


<pt(1—p) (Sn + En t+Vn +18 + 18! + 152 +121 + 1h2 + 1R8 + Rn) 


Finally, we can say that 


s Ss h h h 
(Soa, Enis Vntis tds Tal tne 1 Fost} EQ. 
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2.2 Vaccination reproduction number R, 


The basic reproduction number, R,, is a fundamental epidemiological param- 
eter that quantifies the contagiousness of an infectious disease. It represents 
the average number of secondary infections generated by one infected person 
in a susceptible population. If R, is greater than 1, then the disease can 
sustain transmission; if it is less than 1, then it is likely to die out. 

In this subsection, we describe the method used to estimate the vaccina- 
tion reproduction number R, (also called a control reproduction number) for 
our SEIR-V model. This method, proposed by Diekmann, Heesterbeek, and 
Metz [10], and further elaborated by Van Den Driessche and Watmough [22], 
involves using the next-generation matrix to determine R, for an ordinary 
differential equation compartmental model. 

Starting from the system (2), we can obtain the nontrivial DFE point as 


follows (is calculated in the next section): 


* mS” yom S™ ) 
Eo =|S ,0, —_——_, 0, 0, 0, 0, 0, 0, 
: ( (72 +0 + H) (0 + 1) (n2 + 70 + #) 


with 
(8 + pt) (N2 + 0 + LB) 
(0+ pw) (M2 +70 +H) +m (8+ LW) +10 


Let X = (E, V, 1%, 191, 782, J, Tho, Ths)? Then the model can be written 
as 


= 


X(t+1)— X(t) 
= F(X (t)) — V(X (t)) 
aS (t) (1% (t) + 1° (t) + 1°? (t)) — (81 + B2 + Ba +3 +p) E(t) 
mS (t) + n3B (t) + nal® (t) — (n2 + y0 +) V(t) 
BiE (t) — (na +71 + Bw) I* (t) 
BoB (t)— (Ar +e +H) 1°? (t) 
Ai T°) (t) — (y2 +) T°? (t) 
B3E (t) — (Ae +p) I" (t) 
A2I™ (t) + ell (t) — (As +78 + 61 + w) I" (#) 
Ast”? (t) — (ya + 62 + w) 1” (t) 


We represent F(X) as the rate of appearance of new infections in each 
compartment, while V (X) signifies the disparity between the rate of individ- 
uals leaving each compartment and the rate of individuals entering the same 
compartment through all other means. It is assumed that each function is 
continuously differentiable at least twice in each variable. So, we have 
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aS(I% + [81 4 182) 
0 


o oc oc oO 2 & 


and 
(81+ 62+ 83+73+ uw) E 
(n2 +70 +H) V — mS — n3£ — nal* 
tna yy ee) — Bee 
Artgetyp)l® — BE 
(yo# pa) 1% — Al 
(A2 + w) I — B3E 
(As +93 + 61 + pw) 1% — Ag” — pI 
(ya + da + pe) T% — Agi™ 
The vaccination reproduction number, Ry, is calculated by next generation 


technique. The F and V matrices at DFE is given as follows: 


00111000 
00000000 
00000000 
00000000 
Eo 00000000 
00000000 
00000000 
00000000 


and 


where 
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Vii = (61+ 62+ 63+73+ 4), | V54=—A1, 


V2.1 = —13; V5.5 = (V2+ 4), 

V31 = —fi, Veo = (A2+ pH), 

Vai = —Pa, V7.4 = —Y, 

Voi = —Ps, V7.6 = —Aa2, 

V2.2 = (m2 +7%0+ BH), Vr7= (As +73 + 61 +H), 
V2,3 = —4, V3.7 = —As, 

V3.3 = (m+ +4); Vas = (ya + 62+ y), 
Vaa=OAityrt+y), 0 otherwise. 


Therefore, the next generation matrix is 


a Py eee | . 
Cee a me ee 
(FV- )is= (mtn+H) +74) ¢ 
a 
Oe te = =) oe 
(FV) 


1.5 2+ (y2 +B)’ 
0 ee 


According to Van Den Driessche and Watmough [22], the vaccination repro- 
duction number, R,, is the spectral radius of the next generation matrix, 


which is given by 


> 


R — gge Bz ita + 11 +H) (12 + H+) + Bi (Ar +e + pr) (yo + fH) + A182 (na +71 + 1) 
° that hrm tim tn te) ter aCr ta) 


where 
St = (0 + 1) (n2 + yo + p) 
(0+ pw) (m2 +70 +u+m) +770" 


2.3 Equilibrium points 


The equilibrium points for the model (2) can be obtained by solving the 


following equations: 


S(t+1)—S(@)=E(t+)—-E()=V(t+)—-Vit)=0 
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I(t +1) — I(t) =I! (t +1) —I™ (t) = 1? (t +1) — I? (t) =0 
MmGg+1)-M™ Os ¢4)-™®O)s=F*¢+)-I ®)=0 
R(t+1)-R()= 


The system (2) has two equilibrium points: the DFE and the endemic equi- 
librium. For strictly positive parameters, the DFE is given by 


S* S* 
Eo = (so. 9.0, 0,0, 0,0, acu ) ? (4) 
(j2 + Yo + 1) (0+ ) (n2 +o + 1) 
with 
_ (0+ 1) (m +70 + #) 
(0+ w) (m2 +7 +6+m) +m0' 


and the endemic equilibrium, 


Ei = (ye ie te a a oe ae (5) 


exists if R, > 1 with 


ge = 1 Kk, K3K4K5 
a B2K3 (Ks + A1) + Bi Ka ks’ 
pe — _ KeKsKiksKeKrKsKo (, 1 
k ee 
vere 1m Kk, K3K4Ks5 
a Ko B2K3(K5 +1) + 6 Kak’ 
i 4 481 \ K3K4K5KeK7Ks3k9 1 1 
Le | 13 K; K R, ’ 
pot = 1,3, K2KaKoKoKr Keo () 1 
aie K R,)’ 
. KyK3KsKeK7KeKo 1 
~. = 1 
p82 K R, ’ 
= KyK3KeKrKaKo 1 
[2 = 1 
prA1 Bo K Page 
m KyK3K,KsK7KeKo 1 
ys = 1 
13 K R, ’ 
PY =p A2B3 , PB2\ KoKsKaksKeKaks (,_ 1 . 
Ke Ky K Ry 
net A283. pB2\ Kok3K4K53KeKo 1 
["s = LAs t 1 
Ke K4 Kk Ry)’ 
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pe — 1 20m Kk, K3K4Ks (2 (n mh) 181 x 2.1 Bo 
a Koko Bo K3 (Ks +1) + 6, K4ks Ko K3 K3 KsKk4 


473 A2B3 4 Ba 4 yarA3s (A283 Pe 
K7\ Ke | Ka] K3K7\ Ke ~ Ka 


, KoKsKaKsKoKrKs ( ae ) 


K Ry 
and 

Ko =(m+z); Ks =(72+4), 
Ky = (61+ 62+ 63t+n3t+p), | Ke=(A2+n), 
Ko =(n2+7+ 4); Ky = (Az +493 + 614+ 4), (6) 
Kk3=(m+u+4), Kg = (ya + 62+ 2), 
Ka=(1+¢+h), Kg =(0+1p), 

K = po k3hyKsKoK7 Kg ho + pnsk3Kh4Ks Ke k7 Kg ko (7) 


t+unsyoK3K4aks Ke k7 Ks + wysb2A39KoK3 Ks Ke 

+u74A2A3 3 Ko K3Kaks + W603 KoK3Kh4Kk5K7 Kg ko 
+3283 Ko hK3h4Ks Ks + uBokoK3K5Keh7 Kg kg 
+1 82KoK3KeK7 Kg ko + wbiyKsksKeKk7Kgko 


+2412 KoK3 Keo K7Ks + Bin Ksaks Ke k7Kgk9 

+1 Pim Kaks Ke K7 Ks + u6ryombKsks Ke k7 Ks 
+p73BopKoK3ks Kok + whi (ma + p) Ka ks Ke k7Kgko 

+8171 (Yo + H) Kak Ke K7 Kg + A283 (01 + uw) Kok 3 Ka Ks Kg ko 
+A2A393 (62 + 1) KoK3k4ks Ko + pBo (61 + 1) KoKk3 Ks Ke Kgko 
+A3 982 (62 + 1) Ko k3K3K6Ko. 


3 Stability analysis 


In this section, we study the asymptotic stability of the equilibrium points 
of our system (2) in which the Jacobian matrix is given by 
FS ih Vall sl we ie a ee) Sg) 


No“ni-n In I*n oN ITN 2 i,7=1,10 
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M1 => —Ko —-aQ@ (12 + [st + Is?) 5 
Maa Sole deh + 1p?) 


M31 =, 

Mo.=—-khi, 

M32 = Ns, 

Ma 2 = Ai, 

Ms,2 = Pa, 

M72 = Bs, 

M,3 = Ne, 

Mp4 = Mp5 = M2,6 = 0Sn, 

Mi 4 = M15 = Mig = —aSn, 
where (Kj);9 


M3,3 = —Ka, 
M3,4 = 14; 
Maa = —Ks, 
Mio0,3 = Yo; 
Ms,5 = —Ka, 
Mos = 1, 
Mio4=%1; 
Me,6 = —Ks, 
Ms,5 = 9, 
M77 = —Ke, 


9 are positive constants given in (6). 


3.1 Asymptotic stability of DFE Epo 


Ms,7 = 2, 
Mio = 72; 
Ms,3 = —K7, 
Mo,8 = As; 
Mio,8 = 73; 
Mo,9 = —Kg, 
Mio0,9 = 4; 
Mi,10 = 9, 
Mio,10 = —Ko, 
0 otherwise, 


150 


Theorem 1. The DFE of the system (2) is locally-asymptotically stable in 


Q whenever R, < 1. 


Proof. The system (2) has a unique DFE described in (4). 


The Jacobian matrix at Eo is given by 


J (Eo) = (Mi); 


where 


i,j=1,10? 
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Mi,1 = —Ko, es Mg,7 = 2, 
M31 =; Mas _ Mio,6 = 72; 
M2,2 = —K1, Ms,3 = —K7, 
; Maja = —Ks3, 
M3,2 = 13; Wash Mo,8 = 3, 
M4,2 = Ai, ; Mi0,8 = 73, 
Ms,5 = —Ka, 
Ms,2 = Bo, Me och Mo,9 = —Ks, 
30 a4 ? 
M7,2 = 83, Mio,9 = 74; 
Mio,4 = 71; 
M1,3 = Ne, Mi,10 = 9, 
. —aKgK2 Meo=—Ks, | 0 
M4 = M15 = M16 = a Mio,10 = —Ko, 
Ko Ko Bis "170 85 =P, 
aKg ih? _ 
Mo2,4 = Mo,5 = M26 = , M77 = —-Ke, Be cacneee 
Kg K2 + K9 + 0 otherwise. 


The characteristic polynomial of J (Eo) is given by, 


PIX) = X94 OX Gk tO 4 Oe Cx” 
+ OyX* + OgX? + CgX? +0,X + Cp. 
If Ry < 1, then the coefficients (C;);_9-9 are positives (Details are provided in 


Appendix 6). According to the criteria of Descartes, all eigenvalues of J (Eo) 


are negative reals or complexes with negative real parts. Hence, Eo is locally 


asymptotically stable in 2. 


3.2 Asymptotic stability of the endemic equilibrium FE, 


Theorem 2. The endemic equilibrium of our system (2), E1, is locally- 


asymptotically stable in Q whenever R, > 1. 


Proof. The system (2) has a unique endemic equilibrium FE in Q, given in 5. 


The Jacobian matrix at E) is given by 
J (Fi) = (Mi5); ;-270 


such as 
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a ‘ee be M13 = No, Ms,7 = A2, 
M1 =—(Ko+a(I* +18 +J%°)), 
- a ee M33 = —Ke2, | Mio. = 72, 
Mo, =a(I* +1 + 1%"), 
M34 = "a, Mg.3 = —K7, 
M31 =, 
Maa =—K3, | Mo = Az, 
Mo2=—-Ky, 
M10,3 = Yo; Mio,8 = 73; 
M32 = 13, 
Ms,5 = —Ka, \ Moo = —Ks, 
Ma,2 = fi, 
Mo = 1, Mio0,9 = 4; 
Ms,2 = Bo, u M 6 
10,4 = 1; 1,10 — 9%, 
M72 = 63, 


Me, = —Ks, | Mio,10 = —Ko, 
Mg,5 = 9, 
M77 =—Kée, | 0 otherwise, 


My 4 = M15 = M16 = —-aS™*, 
Mo,4 = M25 = M26 = aS**, 


where 
geal Kk, K3k4ks 
a B2K3 (Ks + A1) + Bi Kaks’ 
7 KW KR RRR 1 
[% = 1 
Trem K R, > 
- Ka hR RRR RS 1 
~. = 1 
[Be K R, > 
si Kok3KeK7Ksk 1 
13" = pr Bo cate (1 x): 


The characteristic polynomial of J (E£)) is given by 
P(X) = X™40y KP 40g RPO X40 X40, XP +O X 4g NEF OLX OX 40H 


If Ry > 1, then the coefficients (C;);_g-g are positives (Details are provided in 
Appendix 6). According to the criteria of Descartes, all eigenvalues of J (£1) 


are negative reals or complexes with negative real parts. Hence, EF; is locally 


asymptotically stable in 2. 


A Sensitivity analysis 


To effectively reduce the transmission and prevalence of COVID-19, it is cru- 
cial to examine and understand the relative significance of various factors 


that contribute to their evolution. The vaccination reproduction number, 
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denoted as R,, is a key parameter that directly affects disease transmission 
control. It measures the average number of cases an infected person gener- 
ates during their infectious period in a population that is fully vaccinated. 
The long-term disease prevalence is represented by the endemic point F. 
Conducting sensitivity analyses on these indicators can help identify the op- 
erational parameters that have a substantial impact on them [5, 20]. This 
knowledge can then be used to develop effective intervention strategies. 

The normalized forward sensitivity index of a variable, f, that is differ- 
entially dependent on a parameter, x, is defined as: 
Of x 


ae da f (8) 


This section focuses on conducting local sensitivity analyses. Initially, sensi- 
tivity indices of R, are calculated, followed by sensitivity indices of FE, using 


parameter values obtained from various sources listed in Table 3. 


4.1 Sensitivity index of R, 


The formula for calculating R, is explicitly given in (3). By using the an- 
alytical formulas provided in reference (8), sensitivity indices for R, can be 
computed and presented in Table 1, using the parameter values from Table 
a 

The sensitivity indices presented in Table 1 are sorted in descending order, 
indicating that the transmission rate a and the recovery rate of infected with 
moderate symptoms yz are the most sensitive parameters. For instance, a 
10% decrease (increase) in the transmission rate a results in a 10% decrease 
(increase) in R,. Similarly, a 10% decrease (increase) in yz leads to a 9.404% 
increase (decrease) in R,. Therefore, it is crucial to accurately estimate 
these parameters, as even a small variation in them can cause significant 
quantitative changes. On the other hand, the lowest sensitive parameters, 
such as the recovery rate of vaccinated individuals yo and the reinfection 
rate of vaccinated individuals 72, have a negligible effect on R,. A 100% 
increase (or decrease) in Yo or 72 leads to only a 0.047043% or 0.041779% 


increase (or decrease) in R,, respectively. 


Tran. J. Numer. Anal. Optim., Vol. 14, No. 1, 2024, pp 136-171 


Keddar and Belhamiti 154 


Table 1: Sensitivity indices of R,. The parameters are ordered from most sensitive to 


least 


Parameters Sensitivity index 
a 1.0000 
V2 -0.9404 
B2 0.6249 
By -0.6125 
V1 -0.0225 
Ay -0.0202 
3 -0.0104 
"4 -0.0040 
Bs -0.0018 
9p -0.0018 
M1 -0.0017 
8 0.00066625 
Yo 0.00047043 
N2 0.00041779 


4.2 Sensitivity index of EF, 


The endemic equilibrium for the COVID-19 system (2), is given by (5). Using 
(8), the same approach can be followed to determine the sensitivity indices 


of each of the endemic equilibrium state variables. 


Table 2 reveals that the most sensitive parameters for all infected states 
are the transition rates (A;);_1 9 3, tates at which exposed individuals become 
infected (3;);1 9.3, and recovery rates (7i);-193,4, respectively. Therefore, 
any variations in these parameters will significantly impact the proportion 
of infectious individuals. The reinfection rate @ is another parameter that 
strongly and positively affects the prevalence of COVID-19, and it is the most 
sensitive parameter for the recovery state. Increasing this parameter of 10 
percent results in a decrease of 9,466 percent in the proportion of recovered 
individuals, leading to an increased risk of COVID-19 spread. Lastly, the 
parameters 74 (vaccination rate of asymptomatically infected individuals) 


and 7, (recovery rate of symptomatic infected) increase the proportion of 
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Table 2: Sensitivity indices of state variables at endemic equilibrium for baseline pa- 
rameter values in Table 3 


” + 


st E* y" 2 1s 8 me me ms Rt 


a al 0.2800 0.2800 0.2800 0.2800 0.2800 0.2800 0.2800 0.2800 0.2734 


By 0.6125 -0.2254 0.6763 0.7746 0.2254 0.2254 -0,2254 0.2254 0.2254 (0.4102 


By -0.6249 0.7206 -0.7206 0.7206 0.2794 -0.7206 0.7206 -0.4531 0.4531 -0.3663 


Bs 0.0018 ——_-0,0069 -0.0069 -0,0069 -0,0069 0.9931 0.9931 0.7256 0.7256 -0.0051 


mh 0 0.0047 —--0.000468 + -0.000468  -0.000468  -0.000468 -0.000468 -0.000468 -0.000468 0.0047 


N 0 -0.0198 -0.2479 0.0020 0.0020 0.0020 0.0020 0.0020 0.0020 _-0,0198 


13 0.0104 0.0045 0.0946 -0.0036 -0.0036 -0.0036 -0.0036 -0.0036 -0.0036 0.0045 


Ny 0.0040 -0.0206 0.7669 0.1495 0.000301 0.000301 0.000301 +~—:0,000301 ~—:0,000301 ~—_—-0.0206 


Yo 0 0.0022 -0.7478 0.0022 0.0022 0.0022 0.0022 0.0022 0.0022 0.0239 


Yn 0.0225 0.0159 -0.7506 0.8341 0.0159 0.0159 0.0159 0.0159 0.0159 0.0369 


Y2 0.9404 0.5741 0.5741 0.5741 -0.4145 -0.4145 0.5741 0.5741 0.5741 0.5800 


¥3 0 0.0036 0.0036 0.0036 0.0036 0.0036 0.0036 -0.7376 -0.7376 (0.0034 


Ya 0 0.0025 0.0025 0.0025 0.0025 0.0025 0.0025 0.0025 0.9775 0.0025 


Ay 0.0202 0.0146 0.0146 0.0146 0.9833 0.0146 0.0146 0.2523 0.2523 0.0146 


A, 0 0.000061 —0,000061 —_-0,000061 1.0001 -0.9998 0.9998 0.000188 0.000188 0.000061 


Ay 0 0.0024 0.0024 0.0024 0.0024 0.0024 0.0024 -0.2445 0.7555 0.0025 


on 0 -0,0058 -0.0058 -0,0058 -0,0058 -0,0058 -0,0058 -0.0177 0.0177 -0.0058 


5, 0 -0.0025 -0.0025 -0,0025 -0.0025 -0.0025 -0.0025 -0,0025 0.0223 -0.0025 


0 0.0535 0.0535 0.0535 0.0535 0.0535 0.0535 0.0535 0.0535 -0.9466 


Q 0.0018 ——_-0,0010 -0.0012 -0,0012 -0,0031 -0,0031 -0,0012 0.2658 0.2658 —_-0,0010 


vaccinated individuals. Therefore, the population becomes immune to the 


disease, ultimately leading to the eradication of COVID-19. 


5 Numerical simulation 


This section aims to provide numerical simulations of the proposed model to 
demonstrate its practicality and effectiveness, and to illustrate our theoretical 
findings. The operational parameters used in the simulations are in Table 3. 

Given that the model is discrete, the transitions between compartments 


are governed by the following rules: 


e For the compartment of exposed individuals, those with a negative 
test result who will be vaccinated are represented by 73, while the 
rest of the exposed individuals are divided into asymptomatic (6; = 
0.55 x (1 — 73 — )), symptomatic (82 = 0.30 x (1 — 3 — w)), and 
hospitalized (63 = 0.15 x (1 — n3 — z)) individuals, according to [18]. 
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Table 3: The values of operational parameters per week 

Parameters Description Range of values Source 
a Transmission rate. 0.002-0.2 [7,14] 
By Asymptomatic infection rate among exposed. 0.4-0.55 [7,14] 
Bo Rate of infection with symptoms among exposed. 0.2-0.3 [7,14] 
Bs Rate of hospitalization among exposed individuals. 0.001-0.15 [7,14] 
vB Vaccination rate of susceptible individuals. 0.00001-0.2 [7,14] 
Nz Co-infection rate of vaccinated individuals. 0.2-0.35 [7,14] 
N3 Vaccination rate of exposed individuals. 0.0001-0.2 [7,14] 
N4 Vaccination rate of asymptomatically infected. 0.1-0.55 [7,14] 
Yo Recovery rate of vaccinated individuals. 0.65-0.75 [7,14] 
V1 Recovery rate of asymptomatic infected. 0.45-0.85 [7,14] 
V2 Recovery rate of symptomatic infected. 0.015-1.0 [7,14] 
¥3 Recovery rate of hospitalized infected in 2"% phase. 0.2-0.75 [7,14] 
Va Recovery rate of hospitalized infected in 3% phase. 0.95-0.98 [7,14] 
Ay Evolution rate symptomatic phase | to 2. 0.6-0.7 [7,14] 
Az Evolution rate hospitalized phase 1 to 2. 0.95-1.0 [7,14] 
A3 Evolution rate hospitalized phase 2 to 3. 0.2-0.8 [7,14] 

Oy Hospitalized COVID-19 mortality, phase 2. 0.01-0.02 (3,7,14] 

b> Hospitalized COVID-19 mortality, phase 3. 0.015-0.05 (3,7,14] 
0 Coinfection rate of recovery individuals 0.55-0.85 [7,14] 

Q The hospitalization rate of symptomatic in phase 1. 0.001-0.35 (3,7,14] 
Lu Natural mortality rate 0.0001-0.002 [7,14] 


The efficacy of the vaccine is represented by yo in the compartment of 


vaccinated individuals, with 72 = (1—-yo — 1) representing the propor- 


tion of vaccinated individuals who are not protected by the vaccine. 


In the compartment of asymptomatic infected individuals, the recovery 


rate is represented by 71, and n4 = (1—7,—2) represents the proportion 


of asymptomatic infected individuals who will be vaccinated. 


For the compartment of symptomatic infected individuals not requiring 


hospitalization, the proportion of individuals entering the second week 


of illness is represented by A1, while y = (1 — Az — yp) represents the 


proportion of individuals whose condition worsens and requires hospi- 


talization. 
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e In the compartment of hospitalized infected individuals in the second 


phase, the recovery rate is represented by 3, A3 represents the propor- 
tion of individuals who remain hospitalized, and the mortality rate due 
to COVID-19 is represented by 6, = 1 — 73 — A3 — pL. 


Finally, in the compartment of hospitalized infected individuals in the 
final phase, the recovery rate is represented by 74, and the mortality 
rate due to COVID-19 is represented by 62 = 1 — y4 — wp. 


5.1 Numerical stability 


In this subsection, we verify the theoretical results obtained above through 


numerical simulations. 


1. The following parameter values (per week) are used: 


a= 0.0025, yo =0.65, y=0.45, 73 =0.2, (9) 
v4=0.95, m=0.2, 3=0.2, 0=0.55, A3=0.78. 


The existence and stability conditions in theorem 1 are satisfied, 
Ry = 0.0015 < 1. 
Therefore, the system (2) has a DFE in Q, 
Eo = (0.6962, 0, 0.1392, 0,0, 0,0, 0,0, 0.1645) . 


Figure 2 indicates under some conditions that the long-term behavior of 
the solution curves of the system approaches the free equilibrium Eo, as 
revealed by numerical simulations. This suggests that the population 


under consideration gradually becomes free of the infection. 


Based on the numerical simulations depicted in Figure 2, it can be ob- 
served that the behavior of the solution curves of the system approaches 
the DFE Ep over a long period of time. This simulation confirms the 


theoretical results outlined in section 3. 


2. The following parameter values (per week) are used: 
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Immunized 


Infected 5 065 


Susceptibles 


Figure 2: The proposed model has only one DFE £po, which is asymptotically stable 
when Ry < 1. 


a = 0.19, yo =0.75, 71=0.85, ¥3=0.75, 74 =0.98, 


(10) 
m = 0.00089, 73 = 0.0089, 6=0.85, A3 = 0.24. 


The existence and stability conditions in Theorem 2 are satisfied, 
R, = 4.5711 > 1. 
There is only one endemic equilibrium point in Q, 


By = (0.21840, 0.03385, 0.00326, 0.018451, 0.015451, 


0.661920, 0.0000507, 0.000068, 0.000017, 0.033088) . 


Numerical simulations in Figure 3 demonstrate that the solution curves 


of the system approach the endemic equilibrium FE, over an extended 
period. 


Our proposed discrete model for coronavirus disease 2019 underwent nu- 
merical processing, which confirmed the existence of two equilibrium points. 


The first is the DFE Eo, which is locally asymptotically stable when R, < 1. 
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Immunized 


Infected Susceptibles 


Figure 3: The proposed model has only one endemic equilibrium £1, which is asymp- 


totically stable when R, > 1. 


The second is the endemic equilibrium F, which is locally asymptotically 
stable when R, > 1. 


5.2 The effect of the most sensitive parameters of the 


model 


It is important to note that the transmission rate a is mainly influenced by 
non-pharmaceutical measures such as physical distancing, the use of masks, 
quarantine and isolation of infected or exposed individuals, as well as limita- 
tions on nonessential gatherings and travel. Medical care is also an important 
factor that can impact the recovery rate yo, ¥1, 3, and ya for COVID-19. 
This includes access to appropriate medical facilities, equipment, and trained 
health-care professionals. The availability of treatments and medications that 
can effectively manage COVID-19 symptoms can also impact the recovery 


rate by reducing the duration and severity of the illness. 


The primary objective of the simulations was to observe how the model 
parameters most affected the transmission dynamics and spread of COVID- 


19. Table 3 values were used in conjunction with varying levels of a, yo, 71, 
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73, and y4, which were increased in stages, ranging from mild to moderate 
and strict rates. The results of the simulations provide insight into the impact 
of control measures on the spread of the disease. By analyzing the effects of 
different levels, policymakers and public health officials can make informed 
decisions regarding the implementation of measures to curb the spread of 
COVID-19. 


Figure 4: Effect of a transmission rate on different states of COVID-19. 


Figures 4 and 5 present the results of numerical simulations conducted 
to investigate the effects of model parameters on the transmission dynamics 
and spread of COVID-19. The simulations were performed using different 
values of a, Yo, 71,73, and ya, and the results provide valuable insights into 


the impact of control measures on the spread of the disease. 


Figure 4 shows that the rate of transmission is a critical factor in deter- 
mining the number of people infected with COVID-19. The graph indicates 
that the higher the transmission rate a, the greater the number of infected 
individuals. When the transmission rate increases from the lowest 0.06 to the 
highest 0.22, the number of infections increases by more than 4 percent, as 
well as the COVID-19 mortality increases by 0.047e — 2 percent. This finding 


highlights the importance of implementing effective measures to control the 
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Figure 5: Effective recovery rates Yo, 71, 73, and ya on the number of COVID-19 infected 


cases. 


spread of the virus, such as social distancing, wearing masks, and limiting 


large gatherings. 


Figure 5 examines the impact of the recovery rates 70,71, 73, and 74 on 
the number of infections. The graph shows that increasing the recovery rates 
results in a reduction in the number of infections. However, the effect of 
Yo, 71,73, and y4 is less pronounced than that of a. When the values of 
Yo, ¥1, 73, and y4 are changed from the highest to the lowest, the difference in 
the number of infections is around 1.5 percent, as well as the mortality from 
COVID-19 decreases by 0.10e — 2 percent. This finding suggests that while 
recovery rates are essential in reducing the spread of the disease, they may not 
be sufficient on their own. A combination of measures, including vaccination, 
testing, and contact tracing, may be necessary to achieve effective control of 


the virus. 


The results presented in Figures 4 and 5 are consistent with the sensi- 
tivity analysis results obtained earlier. The sensitivity analysis revealed that 
the transmission rate a had a significant impact on the number of infections, 


while the recovery rates Yo, 71,73, and y4 had a smaller effect. The simula- 
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tions conducted in Figures 4 and 5 confirm these findings and provide further 
evidence of the importance of controlling the transmission rate. 

Overall, the numerical simulations provide important insights into the 
transmission dynamics and spread of COVID-19. By varying the values of a, 
Yo, 71,73, and ya, the simulations show how different control measures can 
affect the spread of the virus. The results indicate that while recovery rates 
are important in reducing the number of infections, the transmission rate is 
the primary determinant of the spread of the disease. 

In conclusion, the numerical simulations presented in Figures 4 and 5 pro- 
vide valuable insights into the transmission dynamics and spread of COVID- 
19. The simulations confirm the importance of controlling the transmission 
rate and suggest that recovery rates alone may not be sufficient in reducing 
the spread of the virus. The findings of the simulations can inform the de- 
velopment of effective control measures to curb the spread of COVID-19 and 


help mitigate the impact of the pandemic on public health and the economy. 


6 Conclusion 


In this research, we have conducted a comprehensive analysis of the factors 
affecting the transmission dynamics of COVID-19, employing both sensitivity 
analysis and numerical simulations. Our results shed light on the fundamen- 
tal mechanisms governing the spread of the virus, offering critical insights for 
pandemic management and public health policy. 

Sensitivity analysis revealed that the transmission rate a and the recovery 
rate of individuals with moderate symptoms (72) are highly influential pa- 
rameters. Small changes in these values can result in substantial alterations 
in the course of the pandemic. This underscores the imperative of precise 
parameter estimation and the significance of control measures that directly 
impact these variables. 

Furthermore, numerical simulations provided a practical perspective, re- 
inforcing the pivotal role of the transmission rate. The simulations demon- 
strated that variations in a have a profound effect on the number of infections 
and, subsequently, the COVID-19 mortality rate. This underscores the im- 


portance of strict public health interventions such as social distancing, mask- 


Tran. J. Numer. Anal. Optim., Vol. 14, No. 1, 2024, pp 136-171 


163 


A study of global dynamics ... 


wearing, and restrictions on gatherings to control viral transmission. While 
the recovery rates yo, 71, 73, and 4 also play a role in reducing infections, 
their impact is less pronounced compared to a. This suggests that a multi- 
faceted approach, encompassing vaccination, testing, and contact tracing, is 


essential to effectively manage the pandemic. 


In sum, our study emphasizes the critical need for tailored public health 
strategies, underpinned by a thorough understanding of the varying impacts 
of different parameters. These findings are pertinent to real-world policy 
decisions, guiding efforts to mitigate the consequences of the COVID-19 pan- 
demic on both public health and the economy. 


Appendix 


The following calculations were performed using the mathematical software 
(Maxima Version 5.42.1). 


The characteristic polynomial of J(Eo) is given by 


P(X) =X 4 CyX? + OX 4X” + OX 4 Cex” 2 OyX 4 OX? 
A GgN? = OX +O: 


So we have 
min2 < (m+H) (n2 +70 + 1) < Koko, (11) 
auS™ (B81 (Ka + Ks) + Bo (K3 + Ks) + A182) (12) 
. :. 4 i 4 1 
=> as (aksks & + m) + Bo K3Ks5 (<3 + m) + Boks --) 
ns 1 1 1 
< aS" (6, K4Ks5 + 82K 3K5 + A182K3) (= + Ki + ral 


1 1 1 
ky, Kk3k4ksRy {| —+——-4+4— 
< Ny Kg h4its (EtetE) 


< (A, K3k4+ kK, K3Ks5 + kK, K4Ks) Ry, 


and 
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1 1 
as* (By + Ba) = as* (kik Ki Ks + Bo K3Ks5 aK) (13) 


1 1 
< a5" (BKK + Bakke + WA) (a+ FER) 


1 1 
Kk, k3K4ksRy 
Seen aaeh (cate) 


< (kK) K3 + K1K4) Ry. 
We pose 


P, = KoK2— mm, 
Py (Rho ORG IIS 

—aS™ (81 (Ka + Ks) + G2 (Ka +Kg4+A1)), 
P3 = (K).K3 4+ 11K.) — aS* (61 + Be). 


The coefficients of the characteristic polynomial are given by 


Op = Ge Ki Ka Kg KK KoK7KaKo (1 — Ry), 


C1 = Ga KK; Koo (Ka KiKs + Po) 


+ Kk, K3K4Ks5K6K7Kg (Kok 9 + K2Ko + Pi) (1— Ry) 


ze SK Koka KiKoKo (Rok RRR ee i 2), 

Cy = Ke Kk7 Kg ko (Ki K3K4 + K1K3K5 + Kk, K4Ks5 + K3K4Ks) (Ko + Ko) 

+ ky, K3Ky4Ks5Kgkg (Ke K7 + K2k7 + KoK7 + KoKe + KoKe) (1— Ry) 
+ Ky, K3K4K5K6k7 (Ko Ko + Kok 9 + Ko Ks + K2Kg) (1 — Ry) 


an SK KaK Kako Ko (Ke + Ky + Ks) (1 _ Ry) + Ke K7KgP\P2 


+ Ge KaK (Ki K3K6K7Kg + K3KsKoK7Ks + KiKsKoK7Ks) 
+ £KoK3KiKo (KsKoKr + KoKoKs + KoK7Ks + KoK7Ks) 
LU 


Ge KK (Ke K7 + KeKs + K7Kg) Po + Ke K7 Ks ko (Ko + Ko) Po 


+ kK, K3K4K; (Ke k7 + Keo Ks + KgK7) Pi (1— R,) 


+ KgKiKoKoKrKsPi + £KoKoKrKsKoPs, 
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C3 = Ki K3K4Ks5 (Kok (Kg + Ko) + (Ko + K7) Kg Ko + a) (1 — Rv) 


| Ks Ke K7 Ks ko ((K1 + K3 + Ka) (Ko + Ko) + K3K4) + Ke K7KegP3P1 
| Ky, K3Kaks (Kg (K7 + Ke + Ko) + (K7 + Ko) Ko) (Ko + Ka) (1 — Rv) 


+ K3k4K5KeKgkog (Ko + Ka) + ki K3K4K5Ke¢Kk7 (Ko + Ko) (1 _ Rv) 
| K3Ka (K7Ks (Ke Ko + Ks Ko) + Keo K7 (Ks Ko + K5Ks)) (Ko + K2) 
+ Gr koks Ki k5Kokg t+ K3K4K5 Kg + Ki K5K7Ks3 + K3K4Ke6Ks) 


+ ae KoKy K7Kg + KeKks + K7K6) P3 + Ke K7 Ks ko (Ko + K2) P3 
} KK K7 Kg (Ki ke + K3K4+ K3K5) + K3K6Kg (K3 + Ka)) 


, 
S* 

_ 
S* 

+ Fe koKy Kg + K7+ Ke) Po+ Ki k3K4ks5 (Kg + K7 + Ke) Pi (1 — Rv) 

+ (Ko Kz (Ke Kg + Kg Ko 4 Kg Ko) + Ko Kg Ksgkog + KoK6K7 (Kg + Kg)) Po 


Ko Ko (Ke K7 Kg (K3 + Ka + Ks) + Ka K5K7Ksg + K3K4K5 Ke) 


Ko Ko (K5 Ko K7 (Ki + K3 + Ka) + K3K4K5K7 + K3K4Ko6K7) 


(K7Ks + Ke Kg + Ke K7) PiP2 + Kg ko (Ko K6 + KoK7 + Ko K7) P2 


| Ks Keo K7 Ks (Kk3+ Ka + Ki) Pi + K7K3Kaq (K5K6 + K5Kg + Koks) Pi 


Cy = K3Ko9 (K7Ks (Ke + K5 + K4) + Kg (Ks Ke + Ky4Kke+ Ks kK4 ) (Ko + Ko 


| Ka (Kg ko (Ke K7 + K5K7 + Ks Ko) + Ks Ko Kz (Ko + Kg) 
| K3 (KeK7 (K5K3 + Kaks + Kiks) + Ka ks Ks (K7 + Ke) 
| Ky (Kg ko (Ke K7 + K5K7 + Ks Ko) + Ks Ko Kz (Ko + Kg) 


t 5 KaKoKsKr (Ke + Ks + Ky + K3 + 1) + kK, K3KkyKs 


) 
yA 
)( 
) (Ko + Ko) 
(1 


| Ky K3K4kKs ((Kg + K7 + Ke) Kg + (K7 + Ko) Kg + Ke Kz) (1— Ry) 


| Ky Kg Ko (kK) K5Ko6 + K3k5Ko + Kaks5 Ko + K3k.Ks + K3K4K¢) 


| K3K4KsKe (K7Kg + K7Kg + KgKo) + K3KoK7Kgko (Ko + Ko) 


| K3Ko (Ks K6K7 + Kako K7 + K5K4K7 + Ke K5K4) (Ko + K2) 


a 
S* 


K2Ko Ke (Kg + Kz) (Ks + Kg + K3 4+ K1) + Ke K7KgkKoP3 


+ Kg ((K7Kg + Ke Kg + KeK7) + (Kg + Kz + Ke) (Ko + Ko)) Po 


(KG 1G? Be) PP + §KaKsKiKo (hehe he i) 
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+ (K7Kgko + KeKkg kg + Ke K7 Kg + Keg k7Ke) (Ko + Ko) P3 
+ Kok (1 Ks + K3k4+ K1Ks + K3K7 + K1K7 + K3K5) Pi 


+ K7 Kg (K3K4+ kK, Ks + K14Ks 4+ KyKo + Ks Ko + K3Ks5) Pi 


4K Ry BGR + Rake RR KK PP) + $-KaKoP2 


_ 
S* 
+ Kk, K3K4Ks5 (Kg + K7 + Ko + Ke) (Ko + Ka) (1— Ry) 


KoKog (Ke + Kh7+ Ks) (Ks (ky + K3+ K4) + Ps) 


Cs = K3 (Ko + Ko) (K7Kekg + (K7Kg + Ke Kg + K7Ke) (Ks + K4)) 
+ Ks (Ko + Ky + Ko + K3 + K4) (Ko (Keg (K7 + Ke) + K7Ke)) 
+ K4Ko (Ko + Ko + K3) (Ks (K7 + Ke + Ks) + K7 (Ke + Ks)) 

) 


+ K4 (Ko + Ko + K3) (Kg (Kz (Ks + Ke) + Ke Ks) + K7 Ko Ks 


+ K3K4Ks (Ko + Ko) (Kg + K7 + Ke) + (Kg + Kz + Ke) PiPs 
+ ((K7Kg + KeKs + Ke K7) Kg + K7 Kg Ke) Ps + Se KaKoPs 

+ (K7Kg + Ke Kg + Ko (Kg + K7 + Ke) + K7 Ke) (Ko + Ko) P3 
+ Kk, K3K4Ks (Ko + Kg + K7 + Ko + Ko + Ko) (1— Ry) + PoPi 
: SKK (Re (Ri Bo at Re) 4 (4 Ky KD) 

+ ((K7 + Ke) Kg + (Ke + K7 + Kg + Kg) (Ko + Ko) + K7Ke) Po 
+ K3Ko (Ko + Ko) (Kg (K7 + Ke) + (Kg + K7 + Ke) (Ks + Ka)) 


+ ((Kg + Kh7+ Ke) Kg) Po + K3Ko9 (Ko + Ko) (K7K¢ + KsKa4) 
de He (FkoKo fe KsP1) (hie EG RS 
+ Ky (Ko + Ko) (K7Kg Ko + Ke Ks (Ko + Kg + K7)) 


| Ky (Ke + Ko) (Kg ko + K7Ko + K7Ks) (Ke + Ks) 

+ SKK (Kl tie eh Rete) 

i 
S* 

+ Ko Kgk7 Ke (Ko + Ky + Ko + K3 4+ Ky + Ks) 


Kok3K4Ko + KyKsK6Ko (Ko + Ko + K3) 


+ (KeKg + Ke Kz) (Ky + K3 + Ky 4+ Ks) Pi 
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+ K5K6K7 Ks (Ko + Ky + Ko + K3 4+ K4) 
+ K3K4 (Ks + Ke + K7 + Kg) Pi 
+ (Ks (K7 + Ke + Kg) (Ky + K3 4+ Ka)) Pr, 


(Kg (K7 + Ke + K54+ Ka+ K3+ Ki) + K7 (Ko +Ks4+ Ka+ K3+Ki)) Pi 


+ (Ks (Ko + Ki + Ko + K3 + Ka) + K3 (Ko + Ko)) (Ks (K7 + Ke) + K7 Ke) 
+ Kok, (Ko (Kg + K7 + Ke + Ks) + Kg (K7 + Ko + Ks) + K7 (Ko + Ks)) 


4: Fe koks (iG 6 Ry Ky 4+ Ke Rg ie + Be) + RR = R,) 


+ K4 (Ko + Ko + K3) (Ko (Kg + K7 + Ko + Ks) + Kg (K7 + Ke + Ks)) 


+ Ko (Ks (Ko + Ki + Ko + K3 + Ka) + K3 (Ko + Ko)) (Kg + K7 + Ke) 
+ K7 (Koki + Kak3) (Ko + Ks) + Kako Kz (Ko + K2) + Kak3 Ko Ks 
+ (Ko (Kg + K7 + Ko) + KeK7 + (Ko + Kg + K7 + Ke) (Ko + K2)) P3 


+ (K3K4 (Ks + Ke + K7 + Kg + Kg) + Ka ks (Ko + K7)) (Ko + Ke) 


+ Ko ky (Ko (Kg + K7 + Ke + Ks) + Kg (K7 + Ke + Ks) + Ko Ks) 
+ K7Kg kg (Ko + Ki + Ko + K3 + Ka + K5 + Ko) + Ki K2K5 Ko 


Ke (Ko + Ki + Ko + K3 + Ka + Ks) (Ko (Kg + K7) + Kg K7) 


+ (Ke (Ks + Ka+K34+ K1)+ Ks (Ka+K34+ Ki) 4+ Kak3) Pi 
+ K3Ks5 (Ko + K7 + Kg + Ko) (Ko + Ko) + Kg (K7 + Ko) Ps 


(Ko + Kg + K7 + Ke + Ko + Ko) P24+P3Pi, 


C7 = K7KgKo + Ke (Ko + Kg + K7) Ks + (Ko + Kg + K7 + Ko + Ko + Ko) P3 


and 


+ (K7K9 + K7Kg + Ko Ks) (Ko + Ki + Ko+ K3+ Ka+ Ks + Ke) 


+ ((Ke + Ks) (Ko + Kg + K7) + K5 Ko) (Ko + Ki + Ko + K3 + Ka) 


+ K3 (Ko + K2) (Ko + Kg + K7 + Ke + Ks + Ka) 4 ge KaKo + Po 


+ (CK, + K4) (Ko + Keo) + K3K4) (Kg + Kg + K7 + Ko + Ks) 


t (Kg +K7+ Ko +K5+Ka+K3+ Ki)Pi, 


Kg 4+ K7+ Ke) (Ko + Ki 4+ Ko+ K3 + Ka+ Ks) + Ki (Ko + Ko) 


B 


+t Ks (Ko + K1+ Ko + K3 +4 K4) + Ko (Ko + K7 + Kg) + Kg (Ko + K7) 


+ Ka (Ko + Ko + K3) + K3 (Ko + Ko) + KeK7+PitPs, 
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Coke At 4 het Rie py et BK: 


Consequently, if R, < 1 and inequalities (11)-(13) are true, we have 
Cie {0,1,...,9} > 0. 

Using a technique similar to the one previously employed, the coefficients 
of the characteristic polynomial for the Jacobian matrix at the endemic equi- 


librium J (£,) can be shown to be positive if R, > 1. 
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